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A HYBRID SYMBOLIC-NUMERICAL APPROACH 
TO THE CENTER-FOCUS PROBLEM 

ADAM MAHDI\ CLAUDIO PESSOA^ AND JONATHAN D. HAUENSTEIN® 


Abstract. We propose a new hybrid symbolic-numerical approach to the center-focus 
problem. The method allowed us to obtain center conditions for a three-dimensional 
system of differential equations, which was previously not possible using traditional, purely 
symbolic computational techniques. 


1. Introduction 

1.1. Background. Determination of the local stability of an isolated singular point for 
a system of ordinary differential equations (ODEs) is one of the fundamental problems 
encountered across various branches of applied sciences and engineering. For a system 

(1) X = f(x), X G M"', 

where f : M” D A — )• is smooth, and xq is a singularity, i.e. f(xo) = 0, the celebrated 

Hartman-Grobman theorem [12] states that the linearization of (1) or equivalently the 
set of the eigenvalues Ai,...,A,i of the Jacobian matrix ZJf(xo) characterizes the local 
qualitative behavior of the trajectories provided that the eigenvalues have non-zero real 
part, i.e. Re{Xj) / 0. In this case, we say that xq is hyperbolic, otherwise we say that it is 
nonhyperbolic otherwise. To establish the local stability of nonhyperbolic singular points, 
higher order terms have to be taken into account. 

One of the simplest and well-known stability questions is the center-focus (or center) 
problem, originally defined for planar polynomial differential systems, i.e., system (1) when 
n = 2 and f is a system of 2 polynomials in M[x] of some degree m. It consists of obtaining 
conditions on the coefficients of f(x) to distinguish between a local focus (see Fig. 1(a)) or a 
center (see Fig. 1(b)), which has been the subject of intensive research (e.g., [59, 68, 72, 73, 
15, 66, 10, 9, 56, 14, 23]). Although the problem is open in its full generality, it has been 
solved for some important subclasses of planar polynomial vector fields. As an example, 
consider the quadratic system defined by 

ii = u -b -b a2UV -b 03 ^^ 

( 2 ) .22 

V = —u -b OqU -b a^UV -b OqV , 
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Figure 1. An example of a stable focus (a) and a center (b). 

where oi,..., ag G M. The center conditions were established by Dulac [19] and Kapteyn [42]. 
It is well-known (see e.g. [73, 56]) that, for system (2), the Bautin ideal ^ is generated by 
the first three focus quantities of this system [7]. Moreover, the center variety V(.^) C 
of the ideal SS has four irreducible components, namely 

N{gS) = Niluam) U U V(/a) U V(/,o„), 

corresponding to Hamiltonian systems, reversible systems, the Zariski closure of those sys¬ 
tems having three invariant lines, and the Zariski closure of systems having an invariant 
conic and an invariant cubic, respectively. 

The center-focus problem can also be defined for higher dimensional systems [8] and have 
recently been studied for a number of three-dimensional families [21, 11, 27, 49, 50, 51]. We 
continue this study here by applying our new symbolic-numerical approach to a three- 
dimensional system presented in Sec. 1.3 with results presented in Theorems 1 and 3. 

1.2. Computational challenges and the new approach. The process of solving the 
center-focus problem for a specific system of differential equations can be divided into 
three steps [13]. First, the computation of certain number, say p G N, of focus quantities 
(also called Lyapunov quantities), which are polynomials in the parameters of the system. 
Second, finding the common zeros of the polynomial system formed by the focus quantities, 
or more precisely the determination of the irreducible component of the variety of the 
ideal generated by the first p focus quantities. Third, for the system restricted to each 
component, one checks if the necessary conditions for the existence of a center can be 
applied. This typically involves the application of the Darboux theory of integrability or 
reduction to the center manifold. 

Techniques for efficient computation of Lyapunov quantities has been motivated both 
by mathematical and engineering problems. Over the years, a number of algorithms have 
been developed [57, 53, 29, 30, 48, 43, 71]. In this work, we used a method (described 
in [21]) for computing the focus quantities for a system in dimension three, which is based 
on the equivalence of the existence of a center and a local analytic first integral in the 
neighborhood of a singular point (more details are provided in Sec. 2). The advantage of 
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this approach is that it allows to avoid center manifold approximation, which is especially 
important since its power series approximation of analytic or even polynomial systems need 
not converge (e.g., see [1, 60, 51]). 

From the computational point of view, the biggest obstacle in solving the center-focus 
problem for a specihc system is the determination of the irreducible components of the 
variety (i.e., solution set) defined by a certain number of focus quantities. The most common 
approach [2, 32, 24] is the application of computer algebra algorithms for computing the 
primary decomposition of the ideal generated by the focus quantities such as Gianni-Trager- 
Zacharias (GTZ) [31] or Shimoyama-Yokoyama (SY) [58], which have been implemented in 
various symbolic packages (e.g. Singular [35], or Macaulay [34]). The computational 
difficulty related with Grobner basis calculation over the field of characteristic zero was eased 
by implementation of modular arithmetics [70, 20, 55], and successfully used in numerous 
problems [22, 36, 24, 67]. Unfortunately, in practice, the application of algorithms that use 
Grobner bases (also with modular arithmetics) is computationally very heavy and the center 
conditions can only be obtained for specific systems with few parameters. In this paper, we 
replace this particular step and find the common zeros of the polynomial systems formed 
by focus quantities using numerical algebraic geometry techniques (for more details, see 
Sec. 3 and the books [ 6 , 65]). The parallelizablity of numerical algebraic geometry together 
with a regeneration based approach [38, 41] and exactness recovery [4] provides a natural 
alternative to Grobner basis methods. In particular, for the hrst time, we are able to solve 
the center-focus problem for a quadratic, three-dimensional system described next. 

1.3. An application. Consider a third-order differential equation of the form 

(3) u = il + ii + u + f{u,u,u), 

where / = f{u,u,u) £ ]R[u, h, ii] is a polynomial of degree m. Following [49], we can 
equivalently write 

(4) u =—V + h{u,v,w), i) = u + h{u,v,w), w =—w + h{u,v,'w), 

where h{u, v, w) = f{—u + w, v — w,u + w)/2, which we call the standard form of system (3). 
Note that the origin of (4) is a nonhyperbolic singularity at which the associated Jacobian 
has two purely imaginary eigenvalues Ai ^2 = =ti and A 3 = —1. Various dynamic aspects of 
systems of the form (4) have recently been considered, including the center conditions [11, 
18, 21, 50], limit cycle bifurcations [69, 51], Lie symmetries [27], and isochronicity [54]. In 
particular, the center conditions on the local center manifold for system (4), where 

(5) h{u, V, w) = aiu^ + a 2 V^ -|- a^uP' + a^uv -|- a^uw + uqvw, 

were studied in [49]. Although it was possible to compute the first eight focus quantities, 
standard symbolic algorithms (e.g. GTZ and SY) were not able to provide the decomposi¬ 
tion of the Bautin ideal into primes for a general six-parameter system, even over the field of 
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non-zero characteristics. On the other hand, the application of our hybrid approach using 
numerical algebraic geometry to decompose described in this paper, allowed us to obtain 
the center conditions for a general six-parameter system (4). 

Theorem 1. The system (4) with h{u,v,w) as in (5) admits a center on the local center 
manifold if and only if one of the following holds: 

(1) ai = 02 = 04 = 0 

(2) Oi — 02 = O 3 = O 5 = 06 = 0 

(3) oi -|- 02 = 03 = 05 = 06 = 0 

(4) oi -|- 02 = 202 — 0.3 4“ 06 = 03 — 04 — 2 o 5 = 2 o 4 -|- 805 ag = 0 

(5) 2oi — 06 = 2 o 2 +05 = 203 — O 5 -|- 06 = 04 -|- O 5 -|- 06 = 0 

( 6 ) Oi — 02 = 202 + 06 = 04 = 05 -|- 06 = 0 

(7) 2oi -|- 02 = 2 o 2 -|- 06 = 4:03 4“ 506 = 04 = 205 — 06 = 0. 

As an easy conclusion, note that the irreducible components of the center variety (i.e. the 
variety of the Bautin ideal generated by the focus quantities) of system (4) for quadratic h (5) 
are vector subspaces of its six-dimensional parameter space, which was conjectured in [49]. 

1.4. Outline. The rest of the paper is organized as follows. Section 2 summarizes focus 
quantities and their computation. Section 3 summarizes the numerical algebraic geometric 
solving approach along with exactness recovery method used to prove Theorem 3 in Sec¬ 
tion 4. Appendix A presents the Dulac-Kapteyn criterion of quadratic planar systems with 
Appendix B summarizing Darboux theory of integrability. 


2. Focus QUANTITIES COMPUTATION IN 

This section is a review of the method described in [21] (see also [49, 51]) for studying the 
center problem on a center manifold for vector fields in dimension three. Let A :[/—)• 
be a real analytic vector field, such that DX{0) has one non-zero and two purely imaginary 
eigenvalues. By an invertible linear change of coordinates and a possible rescaling of time, 
the system of differential equations u = A(u) can be written in the form 

it = —v + P{u,v,w) 

(6) V = u + Q{u,v,w) 

w = fiw + R{u,v,w), 

where ,0 is a non-zero real number. Let A = {—v + P)d/du + {u + Q)d/dv + {j3w + R)d/dw 
denote the corresponding vector field. A local first integral of system ( 6 ) is a nonconstant 
differentiable function H defined in a neighborhood of the origin in mapping into M that 
is constant on trajectories of ( 6 ), equivalently, H satisfies 

f)zj f)M f)M 

XH := (-„ + P)- + (u + <3)- + (/,. + ii)- . 0 


( 7 ) 
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sufficiently close to the origin. A formal first integral for system ( 6 ) is a non-constant formal 
power series H m u, v and w such that when P, Q, and R are expanded in power series at 
the origin, every coefficient in the formal power series in (7) is zero. 

It is well-known that system ( 6 ) admits a local center manifold at the origin, 

e.g., see [44, Thm. 5.1]. One of the main tools for detecting a center on a center mani¬ 
fold is the following theorem (see, e.g., [ 8 , 21 ]). 

Theorem 2. The following statements are equivalent. 

(a) The origin is a center for X 

(b) System ( 6 ) admits a local analytic first integral at the origin. 

(c) System ( 6 ) admits a formal first integral at the origin. 

In fact, a real analytic local hrst integral from statement (b) (as well as a formal first 
integral from statement (c)) can always be chosen to be of the form H{u, v, w) = ■ ■ 

where the dots mean higher order terms in a neighborhood of the origin in M^. 

The equivalence of statements (a) and (b) is called the Lyapunov Center Theorem with 
a proof presented in, e.g., [ 8 ]. By this theorem, we can restrict our efforts to investigate the 
conditions for the existence of a first integral H which is equivalent to determine necessary 
and sufficient conditions for the existence of a center or a focus on the local center manifold. 

From now on, we assume that P, Q and R in ( 6 ) are polynomials. We begin by introducing 
the complex variable x = u iv. The first two equations in ( 6 ) are equivalent to a single 
equation x = ix + ■ ■ ■, where the dots represent a sum of homogeneous polynomials of 
degrees between 2 and n. Let x denote the complex conjugate of x. We add to this equation 
its complex conjugate, replacing x everywhere by y which is regarded as an independent 
complex variable and replacing w hy z simply as a notational convenience. This yields the 
following complexification of ( 6 ): 

n 

X = ix + OpqrX^y’^z'^, 

p+q+r=2 

n 

(8) y = -iy+ ^ hp^rX^y^^z'', 

p-\-q-\-r=2 

n 

Z = fiz+ ^ CpqrXPy’^Z^, 

p+q+r=2 

where bgpr = dpqr and Cpqr are such that Y^^+q+r =2 CpqrX^x'^w^ is real for all x G C and 
w £ M.. Let X be the corresponding vector held of system ( 8 ) on C^. Existence of a hrst 
integral H{u,v,w) = + v'^ + ■ ■ ■ for system ( 6 ) is equivalent to the existence of a hrst 

integral for system ( 8 ), denoted again by H, of the form 

H{x,y,z) = xy+ ^ Vjkix^y^z^. 

j+k+£=3 


( 9 ) 
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We now investigate the existence of a first integral H for system (8) by computing the 
coefficients ol XH and equating them to zero. When H has the form (9), the coefficient 
Qjki of in XH can be calculated explicitly (see [21]). Except when j = k and ^ = 0, 

the equation gjki = 0 can be solved uniquely for in terms of the known quantities 
with a + f3 + 'y<j + k + i. A formal hrst integral H thus exists if qkko = 0 for all K G N. 
An obstruction to the existence of the formal series H occurs when the coefficient qkko is 
non-zero. This coefficient is the focus quantity and it can be expressed as 

2K-1 2K-2 

(10) 9KK0 = E {j UK-j+l,K-k,0 + kbK-j,K-k+l,o) Vj^k,0 + E (^K-j,K-k,0 Uj^k.l, 

j+k—2 j-\-k—2 

j>0,fc>0 j>0,fc>0 

where we have made the natural assignments uno = 1 and = 0 for a + P + 'y = 2 but 
{a,j3,j) (1,1,0). We know ^rno = 0 and <7220 is uniquely determined, but the remaining 

ones depend on the choices made for vkko^ K G N, > 2. Once such an assignment is 
made, H is determined and satishes 

XH{x,y,z) = g 22 Q{xyf + g^ 2 ,o{xyf ^ -• 

It is known that if at least one focus quantity is non-zero for a choice of vkko^ then the 
same is true for every other choice of the vkko- The vanishing of all focus quantities, i.e., 

(11) 9KK0 = 0 for K >2 

is both a necessary and sufficient condition for the existence of a center on the center 
manifold, otherwise there is a focus (see [21]). 

By Hilbert’s basis theorem, there exists Kq > 2 such that the set of solutions of gxKO = 0 
for all 2 < AT < Kq is equivalent set defined by an inhnite system (11). Since such a Kq is not 
known a priori, we will apply an iterative approach that solves gKKO = 0 for 2 < K < M + 1 
given the solution set of gKKO = 0 for 2 < AT < M. Without knowing Kq, solving using 
any M >2 does always yield necessary conditions. 

3. Numerical algebraic geometry 

Symbolic methods, such as Grobner basis techniques, take an algebraic viewpoint for 
solving systems of polynomial equations. In broad terms, they manipulate equations to ob¬ 
tain new relations describing the solution set. An alternative approach is to use a geometric 
viewpoint which manipulates solution sets. Following a numerical algebraic geometry ap¬ 
proach, solution sets are represented by witness sets that we discuss below. A more detailed 
comparison of symbolic and numerical approaches is provided in [3]. 

The field of numerical algebraic geometry grew out of the use of homotopy continuation 
for computing isolated solutions to a system of polynomial equations. We will first briefly 
explain using basic homotopy continuation on a polynomial system F : —)> C^, that is, 

F{x) = 0 dehnes a system of N polynomial equations in N variables. The idea is to select 
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another polynomial system G : —?• related to F such that G{x) = 0 is “easy” is to 

solve. The simplest example is Gi{x) = x^' — 1 where di = deg Fi, but there is a wide range 
of constructing the so-called start systems which exploit structure in F (see [ 6 , 65] for a 
broad overview). Let S C denote the set of isolated nonsingular solutions of G = 0. 

The next step is to construct a homotopy H : x C —?■ connecting G and F, say 

H{x, t) = F{x) • (1 — t) + 7 • t • G{x) 

for a randomly selected 7 G C. For each s G S, the homotopy H defines a solution 
path Xs(t) such that Xs(l) = s and H{xs{t),t) = 0 with the goal of computing the endpoint 
Xs(0) = limj^o+ Xs{t). In fact, this limit is either a point in C'^ which must be a solution 
of F = 0 or the path is said to be diverging to infinity. By differentiating H{xs{t),t) with 
respect to t, one obtains the Davidenko differential equation 

JxH{Xs{t),t) ■ Xsit) = -JtH{Xs{t),t). 

By including the randomly selected 7 , called the “gamma trick,” the Jacobian matrix J^H 
is invertible along the path for t G ( 0 , 1 ] with probability one and thus one can use predictor- 
correct techniques to track the solution path Xs{t) starting at Xs(l) = s in order to approx¬ 
imate Xs(0). We refer the interested reader to [ 6 , 65] for more details about path tracking 
and using endgames to estimate Xs(0). In the end, the set E C of convergent endpoints 
of all the paths Xs{t) for s G 5* is a superset of the isolated nonsingular solutions of F = 0. 

We now turn our attention to computing the solution set of F = 0, denoted V(F) C C^, 
for a polynomial system F : —?• C”. Geometrically, V(F) can be decomposed into a 

union of irreducible components V(F) = This corresponds algebraically to a prime 

decomposition of the radical ideal generated by F, namely ^^I{F) = Numerical 

algebraic geometry describes an irreducible decomposition of V(F) by computing a witness 
set for each V), called a numerical irreducible decomposition. 

Suppose that V is an irreducible component of V(F) for some polynomial system F. 
A witness set for V is the triple {F, £,IF} where C C is general linear subspace of 
codimension d = dimF and W = V (1 V(£) so that jIFj = degF. Here, the definition 
of general means that £ intersects V transversely, which is a Zariski open condition on 
the Grassmannian of codimension d linear subspaces in C^. The books [ 6 , 65] provide 
for more information about witness sets including performing computations on irreducible 
components which have multiplicity > 1 with respect to F. 

A witness set for an irreducible component V C facilitates additional computations 
that can be performed on F. Of particular interest to the problems discussed in this article 
include the recovery of exact polynomials that vanish on F, determining the existence of 
real points in F, and intersecting F with another solution set. 

With an input polynomial system with exact coefficients, e.g., in Q, one often would like 
exact output. Although the internal computations and witness sets rely upon numerical 
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approximations, there exist techniques for recovering exact answers which can then be 
verihed using exact symbolic methods, which is typically computationally inexpensive. For 
the problems at hand here, we use the exactness recovery technique described in [4] which 
uses a sufficiently accurate numerical approximation of a sufficiently general point on V to 
compute polynomials with integer coefficients that vanish on V. This method is based on 
using a lattice-base reduction technique such as LLL [45] or PSLQ [25]. 

In many applications, only real solutions or components which contain real points are 
of interest, which is the case here. The approach of [37] uses critical points conditions of 
the distance function to determine if V, represented by a witness set, contains real points. 
If y n = 0, then we can disregard this component from further computations. 

We conclude this section by describing the intersection approach built from witness sets 
which is used in the subsequent section. For this situation, we consider a sequence of 
polynomial systems of interest, namely = {fi, ■ ■ ■, fk} for k > 1. Given witness sets for 
the irreducible components of V(Ffc), our goal is to compute witness sets for the irreducible 
components of V(Ffc_|_i) = V(Ffc) n V(/fc+i). For consistency, we assume V(Ffc) C since 
one can easily adjust the methods to work on projective space which will arise below since 
each qkko is homogeneous in ai,... ,06, i.e., V{gKKo) is naturally a hypersurface in P®. 

For the base case, we need to decompose the hypersurface V(Fi) = V(/i), which can be 
readily performed, e.g., via [64]. 

Now, suppose that we are given witness sets for the irreducible components I4q,..., Vk^n^. 
of V(Ffc). For each j e {1,... , Uk}, we need to compute I4j n V(/a:+i) using the provided 
witness set for I4j, say {F/^, Ckj,Wkj}- Clearly, if vanishes identically on Vkj, we 
know that 14j is an irreducible component of V(i4_|_i). Thus, we shall assume that 14,j 
is not contained in the hypersurface V(/a;+i) so that 14 j n V{fk+i) is either empty or 
consists of irreducible components of dimension one less than 14,j- With this assumption, 
if d := dim 14,j is zero, we know 14,j Pi V(/fc+i) = 0. Thus, we also assume that d> 0. 

We compute Vkjr\V{fk+i) using a regenerative intersection approach developed in [40, 41] 
which builds on the diagonal intersection [63] and the regenerative cascade [39, 38]. It 
can be performed using Bertini [5]. To perform this computation, we select a general 
hyperplane Ti and codimension d — 1 linear space K, so that Ck^j = % r\ K,. Our first 
goal is to compute the finite set of points 14 ,^ n /C n V{fk+i) given II4,j = 14 ,j n /C n 
li g = degfk+i, we select general hyperplanes T-Li,... ,'Hg and compute Vkj H JCnTde for 
i = 1,..., g hy standard homotopy continuation from YkjCifCnT-l. Thus, we have computed 

Vk,j n /C n 

which, again by standard homotopy continuation, can be used to compute I4,jn/CnV(/fc+i). 

Now, to compute witness sets for the irreducible components of 14,^0V(/fc+i), which will 
have the form {Fk+i,JC, •}, we simply need to partition the set of points I4,j CiJCn V{fk+i) 
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into subsets corresponding to distinct irreducible components. This is accomplished by 
using random monodromy loops [ 61 ] with the decomposition certified by the trace test [ 62 ]. 
In short, the idea is to move the linear space /C in a random loop in the Grassmannian and 
observe which points are connected by such paths: the points in I 4 j n/Cn V(/fc+i) which are 
path connected by such random loops must lie on the same irreducible component. Thus, 
monodromy loops yield necessary conditions. The (linear) trace test yields a sufficient 
condition which follows from the fact that, as the linear space 1 C is moved in parallel, the 
centroid of points arising from a union of irreducible components must move linearly. One 
can read about further details of computing such decompositions in [6, 65 ]. 

4 . Center conditions for a three dimensional quadratic system 

The following provides a proof of Theorem 1 . We note that, without loss of generality, 
we can always assume that either ag = 0 or og = 1 . The latter follows immediately by the 
change of variables {u,v,w) i-A (x/ag, y/ag, z/ag) and rescalling of time dt = ogdr. Thus, 
the seven cases in Theorem 1 can be split into ten cases, five each for ag = 0 and ag = 1 . 
After showing these ten cases, we then describe how the seven cases of Theorem 1 follow. 

Theorem 3 . Consider system ( 4 ) with h(u,v,w) as in ( 5 ). 

The system ( 4 ) with ag = 0 admits a center on the local center manifold if and only if 
one of the following holds: 

(a) ai — 02 = 03 = as = 0 

(b) Oi + 02 = 03 = Os = 0 

(c) Ol = 02 = 04 = 0 / 

(d) Ol + 02 = 2ai + 03 = 601 — 04 = 4 ai + 05 = 0 

(e) Ol = 02 + 03 = 202 — 04 = 202 + Os = 0. 

The system ( 4 ) with og = 1 admits a center on the local center manifold if and only if 
one of the following holds: 

(f) Ol = 02 = 04 = 0 

(g) 2 ai - 1 = 04 + Os + 1 = 2 o 2 + Os = 203 - os + 1 = 0 

(h) 2 oi T 1 — 2o2 T 1 — 04 = Os T 1 — 0 

(i) Ol + 02 = 402 - Os + 3 = 602 + 04 + 5 = 202 - 03 + 1 = 0 

(j) 4 ai — 1 = 202 + 1 = 4:03 + 5 = 04 = 2 as — 1 = 0 . 

Necessary conditions. 

We first consider og = 0 and take (oi,... ,as) G IP^. Using the notation from Section 3 , 
12(^2) and 12(^3) are irreducible of codimension 1 and 2 of degree 2 and 8, respectively. Now, 
12(^4) has codimension 3 and decomposes into the following components: 5 linear spaces, 3 
of multiplicity 1 and 2 of multiplicity 3 , and an irreducible algebraic set of degree 39 . The 
three linear spaces of multiplicity 1 are (a), (b), and (c). The other two linear spaces are 
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complex conjugates of each other with their union is defined in by 

Cll + 0,2 — 4fl2 “1“ ~ ~ 0- 

Since the real points on this union are contained in (c), we only need to further investigate 
the degree 39 component, denoted ^ 4 ^ 6 , which is not contained in V(g' 55 o). Regenerating 
from Xifi to compute yields 189 distinct points in of which 19 correspond 

to real points. There are 14 real points that do not lie on (a), (b), or (c) of which only 2 
satisfy gmo = 0, namely (d) and (e). We note that (e) has multiplicity 2 with respect to F 5 . 

We next consider ag = 1 and take (ai,..., 05 ) G C^. Similar to the case above, V[F 2 ) and 
12 (^ 3 ) are irreducible of codimension 1 and 2 of degree 2 and 8, respectively. Also, V(F 4 ) 
has codimension 3 and decomposes into the following components: 3 linear spaces, one 
having multiplicity 1, namely (f), with the other 2 having multiplicity 3, and an irreducible 
algebraic set of degree 41. As above, the two linear spaces of multiplicity 3 are complex 
conjugates of each other with their union defined in C® by 

0 \ ~\~ O 2 — 4(^2 T ^4 ~ 2 u 2 T 040 ^ = 2tt2U5 — ^4 ~ ^5 T 1 — 0 . 

Since there are no real points on this union, we only need to further investigate the degree 41 
components, denoted ^ 4 ^ 4 , which is not contained in V(g' 55 o). Regenerating ^ 4^4 yields 4 
irreducible components of ^ 4^4 n V(( 755 o) not contained in (f) or the hyperplane a| + 1 = 0 . 
Three of these are the lines (g), (h), and (i) with the fourth being an irreducible curve of 
degree 244, denoted X^^ 4 , not contained in 12 ( 5660 )• Regenerating A 5_4 yields 71 distinct 
real points not contained in the hyperplane 05 + 1 = 0 nor satisfying (f), (g), (h), or (i). Of 
these, only one satisfies 5770 = 0 , namely (j). 

Sufficient conditions. 

Cases (a) and (b). If the condition (a) (resp. (b)) holds, system (4) reduces to 

u = —V + aiv? + a 2 v‘^ + a^uv, 

V = U + + 02^^ + OiUV, 

w = —w + aiv? + a 2 v‘^ + a^uv, 

with 02 = ai (resp. 02 = —ui). Note that by Theorem 2, it is enough to show that this 
system admits a local analytic first integral at the origin. Since the first two equations are 
decoupled from the third we only need to show that 

(12) u = —V + oiu + 02 V + 04 UV, i) = u + oiu + 02 V + a 4 uv, 

admits a local analytic first integral. In fact, if 04 7 ^ 0 and 02 = 04 , system (12) has the 
inverse integrating factor 

V{u, v) = -04 + 04 (04 + 2 ai) {x - y) + oi (04 + 2 ai)^ (x^ + yx + y^) . 

As 12(0,0) = — 04 , it follows that system (12) has a first integral defined at the origin. If 
04 = 0 and 02 = oi, applying Theorem 4(ii) with o = c = oi, b = d = —oi, A = 2oi and 
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B = —2ai, we have that (12) has a center at the origin and so it is integrable. The case 
02 = —oi (i.e. case (b)) is analogous, since 

V{u^v) = 1 + (2ai — 04 ) X + (2ai + 04 ) y — aia^x'^ — a^^xy + aia^^y^, 

is an inverse integrating factor for system ( 12 ), which is also nonzero at the origin. 

Case (c). In this case system (4) becomes 

u = —V + + a^uw, 

V = u + + a^uw, 

w = —w + + a^uw. 

Note that w = 0 is invariant and is a center manifold for this system. Moreover, the 
restriction of the associated vector field to w = 0 gives rise to a linear center. 

Case (d). For 02 = —oi, 03 = —2ai, 04 = 6 ai and 05 = —4ai the vector field associate to 
system (4) has the invariant algebraic surface F{u, v,w) = w + ai{u — u)^ — 2ai{v — w)"^ = 0 
with cofactor K{u,v,w) = —1. Since F = 0 is tangent to u) = 0 at the origin, it is a 
center manifold for this system. To determine the dynamics on it first we use the change 
of coordinates (tt, v,w) 1 -^ {x + z,y + z, z) that transforms the system into 

X = -y, 

(13) ij = x + 2z, 

z = —z + aix‘^ + Qaixy + ^aixz — aiy‘^ + 4:aiyz. 

The center manifold F = 0 in the new variables is given by F(x, y,z) = z + ai{x — y)‘^ — 
2aiy‘^ = 0. The restriction of system (13) to F = 0 is given by 

X = —y, V = X — 2 aix ‘^ + Aaixy + 2 aiy ‘^. 

Since this system has the following inverse integrating factor (nonzero at the origin) 
V{u,v) = 1 — 4ai {x — y) + Aaf {x^ — 2 xy — y'^) 
thus in this case system (4) has a center on the center manifold. 

Case (e). For ai = 0, 03 = — 02 , 04 = 2 a 2 and 05 = — 2 a 2 system (4) has the invariant 
algebraic surface F(m, v, w) = w — a 2 {y—z)‘^ = 0 with cofactor K{u, v, w) = —1. Since F = 0 
is tangent to tc = 0 at the origin, it is a center manifold for this system. To determine the 
dynamics on it first we use the change of coordinates {u, v, w) i-A {x,y + z, z) that transforms 
the system into 

X = -y - z + a 2 y‘^ + 2 a 2 xy + 2 a 2 yz, 

(14) y = x + z, 

z = -z + a 2 y‘^ + 2 a 2 yz + 2 a 2 xy. 
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The center manifold T = 0 in the new variables writes as F{x, y,z) = z — a 2 y^ = 0. The 
restriction of system (14) to -F = 0 is 

X =-y + 2 a 2 xy + 2 aly^, ij = x + a 2 y‘^. 

This system is invariant by the change of variables (x, y, t) i—(x, —y, —t) so that it has a 
center at the origin. Hence, system (4) restricted to (e) has a center on the center manifold. 

Case (f). In this case system (4) becomes 

ii = —V + + a^uw + vw, 

V = u + + a^uw + vw, 

w = —w + azw‘^ + a^uw + vw. 

It is clear that the plane w = 0 is invariant and is a center manifold for this system. 
Moreover, the restriction of the associated vector field to w = 0 gives rise to a linear center. 


Case (g). If ai = 1 / 2 , 03 = —02 — 1 / 2 , 04 = 2 a 2 — 1 and 05 = — 2 a 2 , then system (4) has 
the invariant algebraic surface F{u, v, w) = —2w + {u — wY + 2a2{v — w)'^ = 0 with cofactor 
K{u,v,w) = —1. Since F = 0 is tangent to in = 0 at the origin, it is a center manifold 
for this system. To determine the dynamics on it first we use the change of coordinates 
(u, v,w) ^ {x + z,y + z, z) that transforms system (4) with conditions (11) into 

X = -y, 

(15) y = x + 22 :, 

i = -^ + x^/2 + (202 - 'i-)xy + a2y^ + 4a2y2;. 

The center manifold F’ = 0 in the new variables writes as F(x, y, z) = —2z + x^ + 2a2y^ = 0. 
The restriction of system (15) to F = 0 is 

9 9 

X = —y, i) = X + X + 2a2y ■ 

As this system is invariant under {x,y,t) 1 —>■ (x,—y, —t), it follows that it has a center at 
the origin, i.e. system (4) under the conditions (g) has a center on the center manifold. 


Case (h). If oi = —1/2, 02 = —1/2, 04 = 0 and 05 = —1. Then the vector field associate 
to system (4) has the invariant algebraic surface F{u, v,w) = w + [(u + rc)^ + {v — w)^ ] /2 — 
(1 + 03 ) = 0 with the cofactor K{u, v,w) = —1 — 2u + 2azw. Since F = 0 is tangent to 
rc = 0 at the origin, it is a center manifold for this system. To determine the dynamics on 
it first we use the change of coordinates {u, v,w) 1 — {x — z,y + z, z), that transforms system 
(4) with condition (12) into 


X 

y 

z 


-y - 2 z + 2(1 + 03)2^ - x^ - y^. 


-z + (1 + 03)2;^ 


xV2 - y‘^/2. 


( 16 ) 


X 
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The center manifold T = 0 in the new variables is given by 

F{x, y,z) = z-{l + a3)z^ + = 0. 

The restriction of system (16) to F = 0 gives rise to a linear center. 

Case (i). For 02 = —ai, 03 = —2ai + 1, 04 = 6 ai — 5 and 05 = —4ai + 3 system (4) 
admits an invariant algebraic surface F{u,v,w) = w + {ai — l)(u — w)'^ + (1 — 2ai){u — 
w){v — rc) + (1 — ai)(i; — w)'^ = 0 with cofactor K{u,v,w) = —1. Since F = 0 is tangent 
to u) = 0 at the origin, it is a center manifold for this system. The change of coordinates 
{u, v,w) 1 -^ {x + z,y + z, z) transforms system (4) under the conditions (i) into 


X 

( 17 ) y 

z 


-y, 

X + 2z, 

—z + + ( 6 ai — 5)xy + 2(2ai — l)xz — aiy'^ + 4(oi — l)yz. 


Again, in the new variables the center manifold is given by F(x,y,z) = z + (ai 
(1 — 2ai)xy + (1 — ai)y^ = 0 and the restriction of (17) to F = 0 reduces to 

X = -y, V = x + 2{l- ai)x^ + 2 ( 2 ai - l)xy + 2 (ai - l)y^. 


l)x^ + 


This system has the following inverse integrating factor (nonzero at the origin) 

V{u,v) = 1 + 4(1 — ai)x + 2 (2ai — 1) y + 4(ai — l)^x^ 
-4(ai - l)(2ai - l)xy - 4(ai - ify'^. 


Hence system (4) has a center on the center manifold. 

Case (j). For ai = 1/4, 02 = —1/2, 03 = —5/4, 04 = 0 and 05 = 1/2 the vector field 
associated to system (4) admits a polynomial first integral 


H (x, y,z) = x^ + y^ — -x^ — ^x^y + 2x‘^z — ^xz^ — y^x 


+y‘^z - ^yz^ - ^x^z + ^x^z^ + ^y^x^ 


3o1qq o 

+ -^y ^ - yz + xyz + -x^ 


—x'^yz — y‘^xz + 2 yxz^ + 



and so it has a center on the center manifold. 


Proof of Theorem 1. 

Case (1). Follows from Cases (c) and (f) of Theorem 3. 
Case (2). Follows from Case (a) of Theorem 3. 

Case (3). Follows from Case (b) of Theorem 3. 

Case (4)- Follows from Cases (d) and (i) of Theorem 3. 


□ 
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Case (5). Follows from Cases (e) and (g) of Theorem 3. 
Case (6). Follows from Cases (c) and (h) of Theorem 3. 
Case (7). Follows from Cases (c) and (j) of Theorem 3. 


Appendix A. Dulac-Kapteyn criterion 

The following theorem provides a criterion in order to determine when a quadratic pla¬ 
nar polynomial system has a center at the origin. It was first proven by Dulac [19] and 
Kapteyn [42], but we present the version given in [16]. 

Theorem 4 (Quadratic Center). The system 

u = —V — bv? — {B + 2c)uv — dv'^, 

V = u + au^ + {A + 2 b)uv + cv'^, 

has a center at the origin if and only if at least one of the following three hold: 

(i) a + c = b + d; 

(ii) A{a + c) = B{b + d) and aA^ — (36 + A)A^B + (3c -|- B)AB‘^ — dB^ = 0; 

(hi) A + bb + f)d = B + ha + be = ac+ bd + 2af + 2df = 0. 

Appendix B. Basic Darboux theory of integrability 

Since, by Poincare theorem, the integrability is closely related to the existence of a center 
on a center manifold (also on the plane), we provide a short overview of the basic notions 
of the Darboux the ory of integrability used in Section 4; for more information see [46, 33] 
and some applications see [28, 47, 52]. 

We say that F = F{x, y, z) G C[x, y, z] is a Darboux polynomial and F = 0 is an invariant 
algebraic surface of the vector field X if and only if there exists a polynomial K{x,y,z) G 
C[x,y,z], the eofaetor of F, such that XF = KF. A the heart of the Darboux theory of 
integrability is the following result [17]: if there exists some number n of pairs {Fj,Kj) for 
which there exists a nontrivial dependency relation ^ OjXj = 0 then is a hrst 

integral of X. 

Consider now the planar system 
(18) x = P{x,y), y = Q{x,y), 

where P,Q G M[x, y], and the associate vector field X = Pd/dx + Qd/dy. Let U be an open 
subset of M^, and let F, D : 1/ —)• M be two analytic functions which are not identically zero 
on U. We say that R is an integrating faetor of this polynomial system on U if one of the 
following three equivalent conditions holds 
dRP _ dRQ 
dx ’ 


dx 


div(i?P, RQ) = 0, XR = -R div(P, Q), 
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where div denotes the divergence. The first integral H associated to the integrating factor 
R can be easily obtained by 

Hix,y) = j R{x,y)P{x,y)dy + h{x), 


where h{x) is chosen such that it satisfies dH/dx = —RQ. Note that dH/dy = RP, so that 
XH = 0. The function V is an inverse integrating factor of the polynomial system (18) on 


U if 
(19) 


8 V ^dV fdP dQ\ 


We note that {V = 0} is formed by orbits of system (18) and R = 1/V defines on 17\{y = 0} 
an integrating factor of (18). We note that if P and Q are quadratic polynomials and the 
origin of system (18) is a center, then there always exits a polynomial function V : —)> M 
of degree 3 or 5 satisfying equation (19), see [26]. 
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